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Abstract 


During  the  High- Resolution  Remote  Sensing  Main  Experiment  (1993),  wave  height  was  estimated 
from  a  moving  catamaran  using  pitch-rate  and  roll-rate  sensors,  a  three-axis  accelerometer,  and  a  capacitive 
wave  wire.  The  wave  spectrum  in  the  frequency  band  ranging  roughly  from  0.08  Hz  to  0.3  Hz  was  verified 
by  independent  buoy  measurements.  To  estimate  the  directional  frequency  spectrum  from  a  wave-wire 
array  the  Data- Adaptive  Spectral  Estimator  is  extended  to  include  the  Doppler-shifting  effects  of  a  moving 
platform.  The  method  is  applied  to  data  obtained  from  a  fixed  platform  during  the  Ris0  Air-Sea  Experiment 
(1994),  and  data  obtained  from  a  moving  platform  during  the  Coastal  Ocean  Processes  Experiment  (1995). 
Both  results  show  that  the  propagation  direction  of  the  peak  wind  waves  compares  well  with  the  measured 
wind  direction.  When  swells  and  local  wind  waves  are  not  aligned,  the  method  can  resolve  the  difference 
of  propagation  directions.  Using  the  fixed  platform  data  a  numerical  test  is  conducted  which  shows  that 
the  method  is  able  to  distinguish  two  wave  systems  propagating  at  the  same  frequency  but  in  two  different 
directions. 
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1  Introduction 


In  the  past  several  years  research  catamarans  were  used  in  various  in  situ  measurements  of  small  scale  air- 
sea  interaction  processes,  including  wind  waves,  atmospheric  turbulence,  subsurface  turbulence,  and  surface 
films  [  Carlson  et  al.,  1988  ;  Bock  and  Frew,  1993  ;  Hara  et  al.,  1994  ;  Bock  et  al.,  1995  ].  Such  platforms  have 
the  advantage  of  having  very  low  profiles,  hence  they  generate  little  disturbance  on  processes  of  interest, 
in  contrast  to  traditional  shipboard  measurements.  Furthermore,  they  do  not  contaminate  surrounding 
water  surfaces  as  long  as  no  oil-based  engines  are  installed.  This  is  an  essential  requirement  for  studies  on 
small-scale  processes  that  are  easily  affected  by  the  presence  of  surfactants.  Examples  of  sensors  mounted 
on  research  catamarans  include  an  acoustic  anemometer,  a  (scanning)  laser  slope  gauge  to  measure  short 
wind  waves,  a  surface  skimmer  to  measure  surface  enrichment,  an  acoustic  current  meter  and  a  hot-film 
anemometer  to  measure  subsurface  currents/turbulence. 

In  addition  to  measuring  various  processes  independently,  some  of  the  recent  efforts  focus  on  the 
interaction  among  different  physical/chemical  processes  at  the  air-sea  interface.  Examples  include  the  effect 
of  gravity  waves  on  gravity-capillary  waves,  and  the  coupling  between  surface  gravity  waves  and  near-surface 
wind  turbulence.  These  studies  require  simultaneous  and  collocated  measurements  of  surface  gravity  waves 
with  other  measurements  of  small  scale  processes.  Although  directional  wind- wave  spectra  alone  have  been 
measured  successfully  from  a  buoy  [  Longuet-Higgins  et  al.,  1963;  Kats  and  Spevak,  1980  ]  or  from  a  fixed 
platform  [Donnelan  et  al.,  1985],  it  is  of  great  interest  to  achieve  similar  measurements  from  a  research 
catamaran  in  order  to  integrate  the  results  with  other  measurements. 

Recently,  Drennan  et  al.  [1994]  have  demonstrated  that  it  is  possible  to  estimate  directional  wind- 
wave  spectra  from  a  swath  ship.  First  they  apply  the  motion  correction  to  the  signals  from  a  wave-wire 
array  mounted  at  the  bow  of  the  ship.  They  then  apply  the  maximum  likelihood  method  to  estimate  the 
directional  spreading  of  wave  energy.  Motivated  by  their  study,  we  have  attempted  to  obtain  directional 
wave  spectra  by  combining  signals  from  pitch-rate  and  roll-rate  sensors,  a  three- axis  accelerometer,  and  an 
array  of  capacitive  wave  wires,  all  mounted  on  a  research  catamaran.  Although  the  basic  principle  of  our 
data  analysis  scheme  is  similar  to  that  in  Drennan  et  al.  [1994],  a  more  careful  motion  correction  scheme  is 
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needed  because  the  motion  of  the  catamaran  is  substantially  more  affected  by  surface  waves.  In  addition, 
we  have  extended  the  Data-Adaptive  Spectral  Estimator  (DASE)  of  Davis  and  Regier  [1977  ]  to  estimate 
the  directional  frequency  spectrum  of  wind  waves  including  the  presence  of  Doppler  shifting.  This  method 
explicitly  makes  use  of  the  dispersion  relation  of  surface  gravity  waves  and  reduces  the  degrees  of  freedom 
compared  with  the  maximum  likelihood  method  used  in  Drennan  et  al.  [1994]. 

In  section  2,  we  describe  the  motion  correction  to  the  wave-wire  signals.  The  scheme  is  applied  to 

s 

signals  obtained  during  High-Resolution  Remote  Sensing  Main  Experiment  (Hi-Res  2)  of  1993  in  section  3. 
The  results  are  confirmed  by  comparison  with  independent  buoy  measurements  of  wind-wave  spectra. 

In  section  4,  we  extend  DASE  to  estimate  the  directional  frequency  spectrum  of  surface  gravity 
waves  including  the  presence  of  Doppler  shifting.  The  scheme  is  first  applied  to  data  obtained  during  the 
Riso  Air-Sea  Experiments  (RASEX)  of  1994  in  section  5.  In  addition,  numerical  experiments  are  conducted 
to  investigate  how  the  method  can  distinguish  secondary  waves  which  propagate  in  a  different  direction 
from  dominant  waves.  Finally  in  section  6,  the  method  of  motion  correction  and  DASE  are  combined 
to  estimate  directional  wave  spectra  from  a  moving  research  catamaran,  using  data  obtained  during  the 
Coastal  Ocean  Processes  Experiments  (CoOP)  of  1995. 
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2  Wave  Height  Measurement  from  a  Moving  Platform 

In  this  section  we  describe  a  method  of  obtaining  a  true  wave  height  signal  by  combining  signals  from  a 
motion  detection  package  and  a  wave  wire.  Both  instruments  are  mounted  on  a  research  catamaran  which 
is  towed  alongside  a  research  vessel. 

2.1  Coordinate  systems 

s 

Two  reference  frames  are  used  in  Section  2.  Since  the  instruments  are  mounted  on  a  moving  catamaran, 
we  define  a  frame  fixed  to  the  platform  as  the  “rotating  frame” ,  with  the  positive  xr  axis  in  the  forward 
direction  and  the  positive  yr  axis  pointing  to  port.  The  zr  axis  is  defined  to  be  along  the  unit  vector 
kr  —  i,.  x  jr>  where  ir  and  jr  sire  the  unit  vectors  along  the  xr  and  yr  axes,  respectively,  i.e.  the  positive 
zr  axis  is  upward.  The  other  coordinate  system  used  for  analyses  is  the  non-rotating  frame  with  axes 
xc>  yc>  and  zc  defined  by  the  long  time  mean  (in  terms  of  wave-induced  motions)  of  the  rotating-frame 
axes.  The  mean  water  level  is  at  zc  =  0.  The  xc,  yCi  and  zc  axes  are  parallel  to  the  unit  vectors  ic ,  jc j 
and  kCl  respectively.  Time  series  for  our  analyses  are  selected  only  if  they  do  not  contain  large  changes  in 
the  speed  and  direction  of  the  platform  and  in  the  mean  current  velocity.  The  non-rotating  frame  is  then 
approximately  equivalent  to  a  fixed  reference  moving  with  a  constant  velocity  relative  to  the  mean  current 
field.  This  criterion  is  verified  by  analyzing  the  velocity  measured  by  the  current  meter.  For  the  time  series 
analyzed,  the  standard  deviations  of  the  xr  and  yr  components  of  the  current-meter  velocity  are  typically 
0.1  m  s”1,  or  about  15%  of  their  means.  We  may  also  characterize  the  non-rotating  frame  by  analyzing  the 
roll,  pitch,  and  yaw  motions  of  the  catamaran  about  the  xc,  yc,  and  zc  axes,  respectively.  The  standard 
deviations  of  these  angles  are  about  4°  for  pitch  and  roll  and  9°  for  yaw. 

2.2  Characterization  of  wave  height 

The  wave  wires  measure  wave  height  relative  to  the  moving  catamaran.  The  catamaran  motion  is  sensed 
by  the  accelerometers  and  angular- rate  sensors.  These  wave  wire  and  motion  package  data  are  combined 
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to  yield  wave  height  relative  to  the  mean  sea  surface,  77,  according  to  the  equation 


T]  —  7]ww  +  Vuiw  —  acc  1} acc  i 

where  rjww  (vertical  component  of  hww)  is  the  sea  surface  height  relative  to  the  base  of  the  wave  wire, 
T)ww~aec  (vertical  component  of  rwti/_acc)  is  the  height  of  the  base  of  the  wave  wire  relative  to  the  ac¬ 
celerometers  (which  is  negative),  and  r)acc  is  the  height  of  the  accelerometers  relative  to  the  mean  sea 
surface  (see  Figure  1).  These  heights  are  expressed  in  the  non-rotating  reference  frame  (xc,yc,zc)  defined 
in  Section  2.L 

Wave  height  is  calculated  using  a  multi-step  procedure  outlined  in  Figure  2.  The  tjww  term  in  (1) 
is  given  directly  by  the  wave  wire.  The  heights  77 WVJ-aCc  and  Tjacc  depend  on  the  pitch  and  roll  angles  6  and 
4>.  These  angles  are  computed  by  integrating  the  pitch  and  roll  rates  measured  by  the  angular-rate  sensors. 
The  acceleration  vector  measured  by  the  accelerometers  is  then  converted  from  the  rotating  frame  to  the 
non- rotating  frame  using  9  and  <f).  The  vertical  acceleration  zc  is  integrated  twice,  resulting  in  the  height 
7]acc,  shown  in  Figure  1.  The  height  T)wm-acc  is  the  vertical  component  of  the  position  vector  rww^acc  from 
the  accelerometers  to  the  base  of  the  wave  wires  (see  Figure  1).  It  is  found  by  transforming  rww-acc  from 
the  rotating  to  the  non-rotating  frame  using  0  and  <j> . 

2.3  Transformation  matrix 

The  accelerometer  signals  and  the  position  vector  fww-acc  are  measured  in  the  rotating  frame  ( xr>yrjzr ). 
These  must  be  expressed  in  the  non-rotating  frame  (xCi  yCy  zc)  to  compute  the  heights  T]acc  and  r}ww-acc * 
A  transformation  matrix  is  used  for  this  conversion.  Consider  a  reference  frame  (x,y,  z)  which  undergoes 
three  separate  rotations  about  each  of  its  own  axes.  In  general,  the  final  orientation  of  the  frame  depends 
upon  the  order  in  which  the  axis  rotations  are  made.  For  the  case  of  rotation  through  angles  6 ,  and  <£, 
in 'that  order,  the  transformation  matrix 
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relates  the  final  orientation  of  the  axes  to  the  original  orientation  (Macomber  and  Fernandez,  1962).  Note 
that  the  angles  0 ,  and  <j>  refer  to  rotation  angles  about  the  z,  y}  and  x  axes  of  the  frame  undergoing 
rotation,  where  positive  angles  are  described  by  the  right-hand  rule.  The  transformation  matrix  for  any 
alternate  order  of  axis  rotations  is  found  through  rearrangement  of  the  matrices  on  the  right-hand  side  of  (2). 
However,  when  rotation  angles  are  small,  the  transformation  matrix  becomes  approximately  independent 
of  the  order  of  rotation.  In  this  case,  (2)  becomes 

1  —  Ip  0 

Tjm  =  $  1  -4>  (3) 

-e  <f>  l 

where  the  angles  <f>,  6,  and  i p  are  now  to  be  expressed  in  radians.  Equation  (3)  is  used  to  convert  from  the 
rotating  frame  (zr,yr,zr)  to  the  non-rotating  frame  {xc,yc,zc).  The  rotating  frame  is  equivalent  to  the 
non-rotating  frame  when  the  latter  is  rotated  through  angles  <p,  9,  and  ip  about  its  axes. 

2.4  Orientation  angles 

The  angular-rate  signals  are  first  converted  from  voltage  to  angular  rates  in  rad  s  1  using  a  linear  calibration. 
The  angular-rate  sensor  signals  contain  “drift”  error  appearing  as  noise  at  low  frequencies.  Thus,  we 
subtract  the  best-fit  linear  trend  from  each  10  minute  time  series.  The  trend  was  on  average  less  than  0.5% 
of  the  standard  deviation  of  the  angles  with  trend  removed.  The  pitch,  roll,  and  yaw  angles  9,  <j>,  and  ip 
are  then  computed  from  the  rate  signals  using  a  numerical  integration  scheme  with  third-order  accuracy. 
The  form  of  this  scheme  for  the  general,  discrete  variable  a,-  is  given  by 

a,-  =  a;_i  +  — -  “U,'  +  a,-_i  +  —  fit— 2  +  0(f iamp)  >  (4) 

Jsamp  L4  ^ 

where  i  is  the  index  of  the  discrete  a,  a  represents  the  time-derivative  of  a,  and  fsamp  is  the  sampling 
frequency.  This  is  the  Taylor  series  approximation  with  third-order  accuracy,  where  the  term  a  is  replaced 
with  the  central  difference  approximation  (see  e.g.,  Hornbeck,  1975  ). 
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2.5  Accelerometer  height 


The  height  of  the  accelerometers  relative  to  the  mean  sea  surface,  r)acc ,  is  estimated  using  the  pitch  and  roll 
angles.  This  calculation  requires  several  steps.  First,  a  linear  calibration  is  applied  to  the  voltage  outputs 
of  the  accelerometers,  yielding  acceleration  signals  xobs  in  m  s-2.  Next,  a  linear  trend  is  removed  from 
each  ten-minute  signal.  The  average  magnitude  of  this  trend  is  1%  or  less  of  the  standard  deviation  of  the 
signal.  The  accelerometer  signals  are  scaled  for  consistency  with  the  magnitude  of  the  gravity  vector.  The 
scale  factor  is  the  ratio  of  9.81  m  s~2  to  the  ten-minute  mean  acceleration  magnitude,  and  averages  0.965 
±  0.004. 

The  vertical  acceleration  in  the  non- rotating  frame,  zacc,  is  computed  using  the  accelerometer 
signals,  xob, ,  and  the  small-angle  transformation  matrix,  according  to 

Zacc  ~ 

Substituting  for  Tjm  using  (3)  yields 

Zacc  —  —9Xob$  +  ‘frijobs  +  *obs  > 

where  iob, ,  yobs,  and  are  the  *P,  yr,  and  zr  components  of  iob,-  The  zacc  signal  is  then  integrated 
twice  in  the  time  domain  using  the  Taylor  series  approximation  of  (4),  yielding  the  accelerometer  height 
i)acc  of  (1).  Prior  to  each  integration,  the  linear  trend  over  ten  minutes  is  removed. 


2.6  Tilt  height 

The  pitch  and  roll  angles  are  also  used  in  the  estimation  of  the  tilt  height  r)ww-acc,  the  height  of  the  base 
of  the  wave  wire  relative  to  the  accelerometer.  This  height  is  related  to  rww-a cc,  the  vector  from  the 
accelerometer  to  the  base  of  the  wave  wire,  through  the  transformation  matrix  T3m: 

f 

T)wu)  —  acc  — 1  — acc)  *  ^  ^ 


where  rww-acc  is  a 


fixed  vector  in  the  rotating  frame.  Equation  (3),  is  used  to  simplify  (7)  to 


=  —0xww—acc  4-  <t>yww—acc  +  zt uw—acc 
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where  xww-acc,  yww-aCc,  and  zww_acc  are  the  xr,  yr>  and  zr  components,  respectively,  of  the  column  vector 


T'tUW—  ccc  • 


2.7  Wave- wire  height 


The  last  unknown  in  the  computation  of  sea-surface  height  given  by  (1)  is  the  wave  wire  height  t)ww.  Its 
calculation  is  analogous  to  the  tilt-height  estimate  of  (7),  and  is  given  by 


Vruw 


(9) 


Here,  hww  is  the  vector  from  the  base  of  the  wave  wire  to  the  intersection  of  the  wire  with  the  sea  surface. 
The  multiplication  is  trivial  since  the  only  non-zero  component  of  hww  is  the  zc-component,  the  wave  wire 
signal  zww.  By  substituting  for  Tjm  using  (3),  (9)  becomes 


3  Application  of  Motion  Correction  Scheme 

3.1  Overview  of  Hi-Res  2  Experiment 

The  main  experiment  of  the  High-Resolution  Remote  Sensing  Program  (Hi-Res  2)  occurred  in  June  of  1993 
and  was  sponsored  by  the  Office  of  Naval  Research  /  Naval  Research  Laboratory.  The  principle  objective 
of  the  program  was  to  investigate  the  accuracy  of  radar  measurements  of  sub-mesoscale  processes  at  the 
ocean  surface.  The  research  catamaran  was  deployed  as  a  means  of  obtaining  in  sttu  measurements  of 
gravity-capillary  waves  (using  a  scanning  laser  slope  gauge),  surface  gravity  waves  (using  a  capacitive  wave 
wire,  and  a  motion  detection  package),  near-surface  atmospheric  turbulence  (using  a  three-axis  acoustic 
anemometer  mounted  at  5  m  mean  altitude)  and  sub-surface  currents  (using  a  three-axis  acoustic  Doppler 
current  meter  at  1  m  mean  depth).  The  catamaran  was  towed  to  the  side  the  R/V  Iselin  at  a  distance  of 
about  20  m  in  order  to  minimize  effects  due  to  the  wind  blockage  and  ship  wake  of  the  vessel.  Data  from 
the  wave  wire  and  the  motion,  detection  package  were  recorded  digitally  at  a  sampling  rate  of  20  Hz. 

The  location  of  the  experiment  was  near  the  coast  of  Cape  Hatteras,  North  Carolina.  During  the 
experiment  the  position  of  the  ship  was  recorded  with  an  on-board  GPS  sensor.  The  two  main  experimental 
sites  are  shown  in  Figure  3  by  large  circles.  One  site  near  36.6’N,  72.8°W  was  called  “off-shore  site”;  the 
other  site  was  called  “near-shore  site” .  Both  sites  were  located  in  the  vicinity  of  the  Gulf  Stream  North  wall 
so  that  frontal  features  could  be  investigated.  In  the  following  section  we  only  show  results  obtained  in  the 
“near-shore  site”.  Also  shown  in  Figure  3  are  the  locations  of  two  discus  buoys  (filled  circles).  These  buoys 
were  deployed  by  a  group  led  by  Dr.  Hans  Graber  of  the  Rosenstiel  School  of  Marine  and  Atmospheric 

Sciences  at  the  University  of  Miami. 

3.2  Results  of  wave  height  frequency  spectra 

Two  examples  of  wave  height  frequency  spectra,  calculated  using  the  method  of  Appendix  A,  are  shown  m 
Figures  4  and  5.  These  spectra  were  computed  using  161  s  time  blocks  and  total  time  senes  lengths  of  121 
and  90  minutes,  respectively.  The  corresponding  ship  tracks  are  indicated  in  Figure  3.  In  Figures  4  and  5 
the  95%  confidence  limits  are  shown  as  vertical  bars.  These  limits  describe  the  statistical  stability  of  the 
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spectral  estimates,  and  are  given  approximately  by 


7pcm  ~  CrmU)  ~  i  -  2 /y/pc,n{f)  ’  (n) 

where  Cnv(f)  and  Cnn{f)  are  the  true  and  estimated  frequency  spectra,  and  P  is  the  number  of  time  blocks 
used  in  the  spectral  calculations.  A  more  exact  form  of  (11)  is  found  by  replacing  P  with  the  number  of 
degrees  of  freedom.  Because  overlapping  time  blocks  are  used  in  the  spectral  calculations,  the  exact  number 
of  degrees  of  freedom  is  unknown  but  is  within  the  range  P/2  to  P. 

Our  frequency  spectra  may  be  compared  to  frequency  spectra  derived  from  the  two  discus  buoys. 
The  spectrum  from  the  buoy  closest  to  the  position  of  the  catamaran  is  shown  as  a  dashed  line  in  Figures 
4,  and  5.  The  frequencies  of  the  swell  peaks  (near  0.1  Hz)  and  the  swell-peak  spectral  densities  agree 
reasonably  well  in  both  cases.  In  Figure  5  local  wind-wave  peaks  near  0.2  Hz  are  indicated  in  both  the 
catamaran  and  buoy  spectra.  However,  the  wind-wave  peak  spectral  levels  are  larger  in  the  buoy  spectra 
than  in  the  catamaran  spectra  by  factors  of  about  2.  This  discrepancy  may  be  partly  attributable  to  the 

differences  in  the  location  and  in  the  time  of  the  two  measurements. 

The  catamaran  spectra  are  affected  by  Doppler  shifting.  A  simple  means  of  accounting  for  Doppler 
shifting  is  to  assume  that  the  wave  field  is  completely  unidirectional  and  moves  towards  the  catamaran 
along  its  forward  axis.  This  assumption  represents  the  largest  possible  Doppler  shift,  from  a  given  intrinsic 
frequency  /  to  the  observed  frequency  /.  The  catamaran  velocity  relative  to  the  mean  current  in  the 
forward  direction  uniquely  determines  the  Doppler  shifting  of  the  waves.  This  velocity  is  calculated  using 
a  long-term  average  of  the  current  meter  velocity  along  the  xc  axis.  The  ordinate  of  Cnn  is  converted  from 
the  observed  frequency  /  to  the  intrinsic  frequency  /  using  the  Doppler-shifted  dispersion  relation: 

/  =  -jT-  (\fi  +  &~Ucbsf/g  -  l)  ,  (12) 

47Tt/o6j  V  ' 

where  g  is  the  gravitational  acceleration  and  U0bs  is  the  speed  of  the  catamaran  relative  to  the  mean 
current.  Furthermore,  the  <%„(/)  values  must  be  scaled  so  that  the  total  wave  energy  remains  the  same 
after  the  transformation  of  (12).  This  scale  factor  is  Sf/Sf,  where  Sf  is  the  frequency  resolution  and  Sf  is 
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the  intrinsic  frequency  resolution,  defined  by 


6f  =  — % —  fx/l  +  87r[/o63(/  +  <5//2)/<7-v'l  +  87r{/()65(/-<5//2)/<7)  . 

4nU0bs  ' 

The  results  of  Figures  4  and  5  (dotted  lines)  show  that  Doppler  shifting  is  quite  small  for  the  frequency 
bandwidth  described. 

3.3  Accuracy  of  wave  height  results 

We  have  found  that  the  region  consisting  of  frequencies  less  than  about  0.08  Hz  is  contaminated  by  noise 
from  the  zr-accelerometer.  Near  0.5  Hz,  there  is  error  due  to  the  resonant  motion  of  the  catamaran  (see 
Hanson  [1996]  for  detailed  analyses).  We  therefore  focus  on  the  region  of  the  frequency  spectrum  from 
roughly  0.08  Hz  to  0.3  Hz,  shown  in  Figures  4  and  5. 

The  voltage  outputs  of  the  accelerometers  and  angular-rate  sensors  are  linearly  proportional  to  the 
quantities  which  they  measure.  They  are  accurate  to  within  3%,  provided  that  the  physical  quantities  do 
not  exceed  the  maximum  measurement  ranges  of  the  instruments.  Time  series  for  which  this  occurred  were 
excluded  from  the  analyses.  The  wave-wires  outputs  are  also  linear  and  accurate  to  less  than  1%. 

The  wave  height  result  is  derived  by  assuming  that  the  catamaran  orientation  angles  are  small,  so 
that  the  transformation  matrix  of  (3)  may  be  used.  The  small-angle  requirement  is  met  for  yaw  angle  by 
selecting  time  series  for  analysis  without  large  catamaran  turns.  The  standard  deviation  of  compass  angle 
was  less  than  10°  for  the  spectra  shown.  The  yaw  angles  had  standard  deviations  similar  to  those  of  the 
compass  angles.  The  pitch  and  roll  angles  were  inspected  after  the  wave  height  calculations.  These  angles 
had  standard  deviations  of  less  than  4°.  The  size  of  the  standard  deviations  of  <j>,  0,  and  i>  indicates  that 
the  error  resulting  from  the  small-angle  form  of  the  transformation  matrix  is  small. 

We  have  examined  the  effect  of  the  transformation  matrix  of  (3),  i.e.  the  tilt  of  the  platform,  on 
the  wave  height  result.  To  investigate  this,  we  have  computed  zr,  the  height  of  the  accelerometer  in  the 
rotating  frame  of  the  catamaran.  This  calculation  is  identical  to  that  for  T]acc ,  with  the  exception  that 
the  step  given  by  (5)  is  omitted.  This  result  is  independent  of  the  angles  <j>  and  0.  We  then  compare 
the  autospectra  of  the  accelerometer  heights  in  the  non-rotating  and  rotating  frames.  In  typical  cases, 
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the  effect  of  the  orientation  angles  increases  with  decreasing  frequency  and  approaches  about  10%  of  the 
spectral  value  at  /  ~  0.1  Hz,  suggesting  that  the  rotational  correction  is  indeed  necessary. 
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4  Estimation  of  Directional  Frequency  Spectrum 

A  wavenumber-frequency  spectrum  S(k,/)  is  calculated  using  the  directional  model  discussed  in  this  sec¬ 
tion.  This  model  generalizes  the  original  Data-Adaptive  Spectral  Estimator  (DASE)  model  of  Davis  and 
Regier  (1977)  to  allow  the  Doppler  shifting  caused  by  a  moving  platform.  The  reference  frame  used  in  the 
following  sections  has  positive  x  to  the  East,  positive  y  to  the  North,  and  positive  z  upwards. 

s. 

4.1  Doppler-shifted  dispersion  relation 

First,  we  describe  the  selection  of  the  discrete  wavenumbers  k  at  which  we  compute  the  wavenumber-fre¬ 
quency  spectrum.  This  is  performed  with  separate  requirements  on  the  wavenumber  magnitude  k  and  the 
direction  9k.  The  direction  of  k  is  discretized  by  using  a  constant  spacing  89  =  2ir/Nke,  where  Nk$  is  an 
integer.  The  wavenumber  magnitudes  k  are  related  to  the  observed  frequency  /  through  the  dispersion 
relation  for  small-amplitude  waves  with  Doppler  shifting, 

2irf =  \fgiz  tanh(«;//)  —  U0b,K  cos  9rei  .  (13) 

Here  g  is  the  gravitational  acceleration,  U0b,  is  the  velocity  of  the  observer  relative  to  the  mean  current 
field,  and  H  is  the  water  depth.  The  relative  wave-observer  angle  9rei  is  defined  according  to 

9rei  =  6k  —  90bs  , 

where  9k  is  the  direction  of  k  and  9ob;  is  the  heading  of  the  observer  relative  to  the  mean  current  field.  For 
convenience,  we  indicate  the  dispersion  relation  of  (13)  through  the  notation  k  =  k {f,9k). 

The  generalized  DASE  method  requires  solutions  to  the  dispersion  relation  (13).  We  first  consider 
the  special  case  of  deep  water.  In  this  case,  tanh (k // )  ~  1  and  solutions  to  (13)  are  given  by 

«  =  ...2  1  2  —  [(<7  -  A*fUob,  COS  Orel)  ±  y/gTg  -  SnfUob,  cosflre,)]  .  (14) 

2t/, ob}  cos  9rei  L 

for  cos  9rei  £  0.  The  trivial  solution  for  cos  9rei  =  0  is 

k  =  {2  «f)2g-1  .  (15) 
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Equations  (14)  and  (15)  present  5  solutions  to  the  dispersion  relation  (13),  listed  below  in  order  of  increasing 
k  for  fixed  /. 

ka  (-)  solution  of  (14)  for  cos  6ret  <  0.  Approaching  waves. 

kp  Solution  of  (15)  for  cos  9rel  =  0.  Waves  moving  perpendicularly. 

klt  (-)  solution  of  (14)  for  cos  9rel  >  0  and  /  >  0.  Long  overtaking  waves. 

Kst  (+)  solution  of  (14)  for  cos0reJ  >  0  and  /  >  0.  Short  overtaking  waves. 

N. 

kr  (+)  solution  of  (14)  for  cos9ret  >  0  and  /  <  0.  Overrun  waves. 

The  notations  (-)  and  (+)  refer  to  the  choice  of  sign  in  front  of  the  radical  in  (14).  Solutions  for  overrun 
waves  are  found  using  negative  frequencies  /  in  (14).  These  five  solutions  are  shown  in  Figure  6  for  the 
example  Urei  =  0.8  m  s-1. 

The  three  smallest-wavenumber  solutions,  ka,  kp,  and  klt,  are  shown  in  Figure  7  for  U0b,  =  0.8 
m  s-1  and  9obs  =  0°.  Together,  these  solutions  form  closed  curves  in  fc-space  for  small  /.  The  two  largest- 
wavenumber  solutions,  k$t  and  kr,  correspond  to  waves  with  intrinsic  frequencies  much  higher  than  the 
observed  frequency  /.  These  solutions  may  be  ignored  in  the  directional  analysis  if  they  correspond  to 
waves  with  negligible  energy  compared  to  the  other  solutions. 

The  solution  klt  is  defined  only  for  /  and  9rei  satisfying 

f  COS  Orel  ^  f •  , 


where 


/.= 


9 

8-rrUobs 


(16) 


Thus,  for  frequencies  /  above  /»,  there  is  a  range  of  9re[  where  no  solution  kLt  exists.  The  curve  corre¬ 
sponding  to  /  =  0.58  Hz  in  Figure  7  shows  such  a  case.  For  /  >  the  ends  of  the  curve  given  by  klt 
attach  to  the  solution  of  kST-  In  practice,  it  is  prudent  to  choose  the  observer  yelocity  such  that  the  waves 
of  interest,  e.g.  the  dominant  waves,  have  frequencies  /  <  /-- 

When  there  is  no  Doppler  shifting  Uob,  =  0,  the  dispersion  relation  (13)  becomes 

5/ctanh(/cf7)  =  (2tt/)2  .  (17) 
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Solutions  for  k  from  (17)  describe  circles  in  A:-space. 

4.2  Theoretical  background 

Methods  of  estimating  the  propagation  direction  of  waves  generally  utilize  an  approximation  to  the  contin¬ 
uous  integral 

f)  =  Q( C,  f)e-a<dC  ,  (18) 

where  S(k ,  /)  is  the  wavenumber-frequency  spectrum  and  £  is  the  lag  vector  between  two  locations  in  space. 
The  term  Q(C,  /)  Is  the  cross-spectrum  between  two  height  time  series  measured  at  locations  separated 
by  lag  £.  The  wavenumber  k  has  magnitude  k  —  2n/a  and  direction  0*,  where  a  is  wavelength  and  Ok  is 
the  wave  propagation  direction.  The  calculation  of  the  wavenumber-frequency  spectrum  applies  to  each 
observed  frequency  /  independently,  and  thus  the  dependence  of  all  spectra  on  /  is  implied  throughout 
the  following  discussion.  The  angle  Ok  is  measured  in  the  horizontal  plane,  is  zero  along  the  x  axis,  and 
increases  in  the  counterclockwise  direction  when  viewed  from  above. 

The  cross-spectrum  Q(C)  is  continuous  in  lag  <f.  But  in  practice,  we  estimate  this  quantity  for 
discrete  lags  using  the  method  of  Appendix  A,  yielding  Q{Qpq)-  The  discrete  lag  Q?q  is  the  vector  from  wave 
wires  p  to  q ,  i.e. 

Cpq  ~  xp  —  xq  . 

The  cross-spectra  are  grouped  in  the  matrix  Q  with  elements  Qpq  —  Q{Cpq)-  Each  of  the  diagonal  elements 
Qpp  for  p  =  1,2, ...  Nww  is  the  autospectrum  of  wave  height  for  the  pth  wave  wire,  and  corresponds  to 
the  frequency  spectrum  Cqq  discussed  in  Section  3.2.  Here,  Nww  is  the  number  of  wave  wires.  The  best 
estimate  of  the  frequency  spectrum,  Qa>  is  given  by  the  average  of  the  autospectra  for  the  Nww  wave  wires: 

Q«  =  jj-Y,Qpp-  (19) 

p=l 

A  simple  approximation  to  (18)  would  be  to  replace  Q(C)  with  the  estimate  Q  and  to  formulate 
an  estimate  of  the  differential  d( However,  given  the  sparse  sampling  in  £-space,  the  finite  area  elements 
corresponding  to  d<f  are  poorly  defined.  We  instead  use  a  model  in  order  to  approximate  (18).  A  general 
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model  can  be  written  as 


sw^EE*'  (20) 

<7=1  P  =  1 

where  apq  is  an  element  of  a  matrix  a  to  be  determined.  The  normalization  factor  h  ensures  that  the  energy 
of  Qa  is  conserved  in  forming  the  estimate  S(k)  and  is  discussed  later.  Next,  we  relate  the  model  result 
S{k)  to  the  true  spectrum  S{k)  of  (18).  We  first  write  the  true  cross-spectrum  Q(<f)  using  the  inverse  of 

(18), 


Q{Q=  [°°  S(I)eir<dl,  (21) 

J  —  CO 

where  we  now  use  the  wavenumber  1.  To  the  accuracy  of  the  estimate  Qpq,  we  may  replace  Qpq  in  (20) 
with  its  true  value,  Q{(pq),  given  in  (21).  This  substitution  yields 


/OO 

-CO 


S{F)w{kJ)dT, 


(22) 


where 


w(k,  f)  =  J2  £  apq(k)eir<”  .  (23) 

<7=1  P=1 

The  function  \h  w(k, f)  j  acts  as  the  window  function  relating  the  estimated  wavenumber-frequency  spectrum 
to  the  true  spectrum.  This  window  must  approximate  a  delta  function  centered  on  l  =  k. 


4.3  Fundamentals  of  the  generalized  DASE  model 

The  means  of  selecting  the  matrix  a  varies  from  method  to  method.  The  DASE  method  of  Davis  and  Regier 
(1977)  was  chosen  for  two  reasons.  First,  it  is  a  data-adaptive  method.  This  implies  that  a  is  selected  with 
reference  to  the  data  Q  so  that  contributions  to  S{k)  from  l^k  are  minimized.  This  is  in  contrast  to  a 
priori  methods,  which  attempt  to  minimize  the  sidelobe  height  of  w(k,l)  without  regard  to  Q.  Second,  the 
DASE  method  constrains  wave  energy  to  lie  on  the  dispersion  curve  k  =  K(f,0k)-  This  dispersion  relation 
is  a  strong  physical  restriction  on  waves.  Other  methods,  such  as  the  Maximum  Likelihood  method  of 
Capon  (1969  ),  require  no  dispersion  relation  and  thus  allow  undesired  freedom  in  the  calculation  of  S(k). 
Three  constraints  are  applied  in  both  the  original  and  generalized  DASE  method: 
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1.  The  window  u/(£,  /)  is  non-negative. 

2.  Energy  at  a  given  frequency  /  is  distributed  along  the  dispersion  relation  k  =  *(/,  0*)* 

3.  The  estimate  S(k)  of  least  error  is  chosen  from  a  restricted  set  of  solutions. 

Constraint  1  is  imposed  by  requiring  that  the  matrix  a  be  Hermitian,  or  equivalently, 

a(k)  =  jT{k)j*{k)  ,  (24) 

where  the  superscripts  T  and  *  denote  the  matrix  transpose  and  complex  conjugate  operations,  respectively, 
and 


7i  72  •  -  •  JNW 


The  window  w{k,l)  of  (23)  may  then  be  expressed  as 

w(kj)=1(kyT(M)YT(k)  , 


(25) 


(26) 


where  e  is  the  row  vector 


(27) 


Equation  (26)  may  be  recast  as 

w(k,l)  =  |7(*)<*r$f  •  (28) 

This  equation  reveals  that  the  window  w(k,l)  is  non-negative,  i.e.  Constraint  1.  This  result  is  useful  when 
combined  with  (22)  in  that  the  estimate  S(k)  with  least  error  corresponds  to  the  minimum  estimate. 


4.4  Generalized  DASE  window  constraint 

The  application  of  Constraint  2  differs  in  the  original  and  generalized  DASE  methods.  In  the  generalized 
method,  we  require  that  the  volume  under  the  window  w(k,l)  integrated  over  the  region  A k  must  remain 
constant  with  9k  for  fixed  /.  The  area  A k  surrounds  the  point  k  in  £-space,  and  is  defined  by  the  limits  of 
integration  in  the  analytical  form  of  Constraint  2.  This  form  is  given  by 

r6k  +  &9j 2  r\i  _  _ 

/  /  Xw(k,l)  d\d6i  =  1  ,  (29) 

Jek-A8/2  JXo 
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where  A  and  $i  are  the  magnitude  and  direction  of  /,  and  A6  is  a  constant  angular  width  parameter  specified 
later.  The  limits  of  the  radial  integral  are  defined  using  the  dispersion  relation  (13),  according  to 


\o  =  K(f-6f/2A)  } 

A  i  =  K{f  +  6f/2,0t) 

where  8f  is  the  frequency  resolution.  Constraint  2  is  needed  because  the  area  A k  varies  with  6k  for  fixed 
/.  This  variation  must  be  accounted  for  in  the  estimation  of  S(k ),  which  gives  the  spectral  density  per  unit 
frequency  and  per  unit  area  in  Ar-space. 

For  convenience,  we  recast  (29)  as 


1(k)L(k)YT(k)  =  l, 


(31) 


where  L  is  a  matrix  with  elements 

A0/2  * Aj  _  „ 

Lpq(k)  =  /  eu<"Xd\dB,.  (32) 

Jek- A0/2  J  A0 

In  practice,  we  integrate  (32)  numerically  using  the  trapezoidal  rule,  with  40  points  across  Ad  and  20  points 
in  the  radial  direction  A  (see  e.g.,  Hornbeck,  1975). 

The  choice  for  the  angular  width  Ad  is  arbitrary.  We  have  tested  the  sensitivity  of  the  results  to  its 
value.  For  the  case  of  three  wave  wires,  we  have  found  negligible  difference  in  peak  height,  peak  location, 
and  angular  distribution  of  energy  for  values  in  the  range  5°  <  A 6  <  75°. 

In  the  case  without  Doppler  shifting,  the  radial  integration  limits  of  (32)  are  independent  of  9k ,  as 
described  by  (17).  In  this  case,  the  constraint  of  (31)  may. be  compared  to  the  constraint  in  the  original 
DASE  method  of  Davis  and  Regier  (1977).  By  neglecting  the  variation  of  l  •  Cpq  in  the  radial  direction  A, 
we  may  derive  an  approximate  form  of  (32): 

Lpq(k)  =  k^x-4)  eil<"det,  (33) 

^  Jdk-se/ 2 

where 

KO  =  K{f-Sf/2,0k)  ^ 

ki  =  n(f  +  Sf/2,9k) 
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The  multiplication  of  the  constant  factor  (/cf  —  Kq)/2  by  60  gives  the  area  of  A k.  By  multiplying  (33)  by 
2k/(k?  —  Kq),  we  obtain  the  constraint  imposed  in  the  original  DASE  method.  Thus,  the  constraint  of  (31) 
modifies  the  original  DASE  constraint  to  account  for  the  radial  thickness  in  £-space  associated  with  the 
dispersion  k  =  «(/,£*).  This  thickness  is  due  to  the  discrete  nature  of  the  observed  frequencies  /,  and  is 
described  by  the  boundaries  given  in  (30). 

4.5  DASE  maximizatiorTproblem 

Next,  we  write  the  expression  for  §(k)  by  combining  (24)  and  (20),  yielding 

S(k)  =  m*)Q7*t(£)  •  (35) 

The  minimization  of  S(k)  is  then  equivalent  to  maximizing  the  ratio 

u(k)  =  JL.  =  7(fc)L(£)7-TW  (36) 

S(k)  l(k)QTT(k) 

The  task  is  now  to  maximize  v  under  the  condition 

^7Q7*T  “  7^7*^  > 

where  dependence  on  k  is  now  implied.  However,  this  maximization  problem  is  not  tractable.  Instead, 
Constraint  3  is  imposed  by  maximizing  u  subject  to  the  condition: 

;  ^Q7*T  =  L7*t  .  (38) 

The  cross-spectral  matrix  Q  is  square  and  invertible.  Multiplying  (38)  by  Q_1  yields 

* n'T  =  q-1l7*t  -  (39) 

This  is  a  standard  eigenvalue  problem  involving  the  given  matrix  Q-1L.  The  wavenumber-frequency 

r 

spectrum  is  then 

S(k)  =  h/umax{k)  ,  (40) 

where  i/mox  is  the  maximum  eigenvalue. 
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The  final  requirement  we  place  on  §(k)  is  that  the  energy  over  the  area  in  k- space  corresponding 
to  a  particular  /  must  match  the  average  autospectrum  energy  Qa,  i.e. 


Nk$ 

Y,  S(km)dAm  =  Qa  ,  (41) 

mrrl 

where  m  is  the  index  of  the  discrete  k  for  fixed  /,  Qa  is  given  by  (19),  and  dAm  is  the  area  given  by 


dA  = 


X  dX  dd\  , 


(42) 


for  a  particular  k  =  km.  Note  that  the  area  dA  differs  from  the  area  in  the  integral  of  (32)  in  that  66  is  the 
angular  resolution  while  A0  is  a  model  parameter.  Combining  (40)  with  (41)  yields 


h  = 


Qa 

J2m=l  dAm/l/max[km)  ' 


(43) 


The  estimate  §(k)  is  then  obtained  by  substituting  (43)  into  (40). 


4.6  Calculation  of  the  directional  frequency  spectrum 

The  result  S  obtained  in  the  previous  section  is  a  three-dimensional  wavenumber-frequency  spectrum. 
(Note  that  the  /  dependence  has  been  omitted  for  simplicity.)  The  spectrum  S(k,  f)  is  now  converted  to 

the  two-dimensional  wavenumber  spectrum  S'(k)  according  to 

l v 

S' (*)■  =  */£$(*, /j)  ,  (44) 

where  N  is  the  number  of  discrete  frequencies  fj.  However,  the  DASE  method  distributes  energy  so  that 
for  a  given  k  there  is  only  one  frequency  bin  fk  which  contributes  non-zero  energy.  The  relationship  then 
simplifies  to 

S'{k)=6fS{k,fk)  .  (45) 

Now  we  may  compute  the  directional  frequency  spectrum  D{f,  9k).  The  wavenumber  k  is  converted 
to  intrinsic  frequency  /  using  the  dispersion  relation  for  small-amplitude  waves: 

f  =  (2jr)-Vff*tanh (*#)  ,  (46) 
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which  simplifies  in  the  deep-water  case  to 


/=  (2tt)  xy/g> c  . 

The  directional  frequency  spectrum  D(f,  6k)  is  calculated  via  a  conservation-of-energy  constraint  expressed 
by 

KdKS'{k)  =  D{f,6k)df  .  (47) 
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5  Application  of  Extended  DASE  Method 

5.1  Overview  of  RASEX  experiment 

The  Ris0  Air-Sea  Experiment  (RASEX)  was  sponsored  by  the  Office  of  Naval  Research  as  part  of  the 
Marine  Boundary  Layers  Accelerated  Research  Initiative,  and  occurred  in  the  Great  Baelt  off  the  western 
coast  of  Lolland,  Denmark  in  the  Spring  and  Fall  of  1994  (see  Figure  8).  The  data  used  for  directional 
analysis  were  acquired  from  a  48  m  tower  located  in  the  near-shore  waters.  This  tower  was  instrumented 
with  an  extensive  array  of  meteorological  sensors.  The  minimum  and  maximum  fetches  at  this  location 
were  about  2  km  and  70  km,  and  water  depth  was  3  to  4  m.  The  results  discussed  in  Section  5.2  are  from 
an  array  of  three  capacitive  wave  wires  extended  from  the  tower  at  the  water  surface.  These  wires  formed 
the  vertices  of  an  equilateral  triangle  with  sides  of  length  1  m. 

5.2  Directional  frequency  spectra  results 

In  this  discussion,  angles  follow  the  meteorological  convention,  with  zero  northward  and  angle  increasing 
in  the  clockwise  direction.  Wind  directions  are  reported  by  the  angle  from  which  wind  blows,  whereas 
wave  angles  describe  the  direction  towards  which  waves  propagate.  Wave-height  frequency  spectra  from 
the  RASEX  experiment  were  computed  using  102  s  time  blocks  and  total  time  series  lengths  of  30  minutes. 
Two  directional  spectra  D{f,6k)  for  RASEX  are  shown  in  Figures  9  and  10.  In  both  figures,  wind-wave 
peaks  are  shown  to  have  angular  locations  closely  aligned  to  the  wind  direction.  The  result  of  Figure  9 
corresponds  to  wind  angle  9W  =  340°,  wind  speed  Uw  =  12  m  s"1,  and  the  wind- wave  peak  is  located  at 
/  =  0.27  Hz  and  9k  =  145°.  Its  peak  value  is  D  =  2.4  m2  Hz"1  rad"1.  About  5  hours  later  (Figure  10), 
wind  speed  decreases  to  Uw  =  8  m  s-1,  and  the  direction  slightly  shifts  to  9W  =  350°.  Correspondingly, 
the  spectral  peak  moves  to  /  =  0.28  Hz  and  9k  =  160°,  and  its  peak  value  is  reduced  to  D  —  1.6  m2  Hz'1 
'rad-1.  The  difference  between  the  wind  and  peak  wave  directions  is  15°,  and  10°  for  Figures  9  and  10, 
respectively.  Since  the  wind  and  peak-wave  directions  may  differ  slightly  in  the  presence  of  variable  fetch 
(Donelan  et  al.,  1985),  as  in  RASEX,  the  observed  correlations  seem  reasonable. 

The  accuracy  of  the  directional  method  is  undoubtedly  a  function  of  the  ratio  of  wavelength  to 
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wave-wire  spacing.  The  peaks  found  at  0.28  and  0.45  Hz  correspond  to  wavelengths  a  =  2-it/k  of  17  and 
8  m,  respectively.  The  spacing  between  any  two  wires  in  the  RASEX  array  is  1  m.  Thus  the  ratios  of 
peak  wavelength  to  wave-wire  spacing  are  17  and  8,  and  the  directional  method  apparently  performs  well 

throughout  this  range. 


5.3  Numerical  experiment 

A  test  is  conducted  to  investigate  the  performance  of  the  directional  method  in  resolving  a  secondary  wave 
peak  with  a  different  direction  than  that  of  dominant  waves.  We  first  synthesize  wave-wire  time  series 
which  contain  two  different  dominant  wave  systems,  by  combining  two  sets  of  experimental  data  taken  5 
hours  apart.  Let  (77  1,172,  %)  be  the  three  wave-wire  signals  used  in  the  calculation  of  the  result  of  Figure 
9,  and  let  (i?i,  rj'2,  r]'3)  be  the  three  wave-wire  signals  used  for  Figure  10.  Note  that  these  results  have  peak 
directions  about  15®  apart.  We  form  a  new  set  of  wave  wire  time  series  by  combining  these  2  sets  according 


to  the  relation 


f)m  r!m  C 


^Irms 
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(48) 


where  c  is  a  constant  in  the  range  0  <  c  <  1,  and  the  indices  m  =  1,  2,  and  3  correspond  to  the  indices 

n  =  3,  1,  and  2,  respectively.  We  define  the  index  of  77  so  that  it  increases  in  the  clockwise  direction.  Thus, 

we  have  shifted  the  wave-wire  data  of  the  17'  set  by  120°  in  the  clockwise  direction.  The  root-mean-square 

heights  T]rrns  and  T)'rms  are  defined  according  to 

»  1/2 

where  N  is  the  number  of  discrete  times  ij  used  in  the  directional  estimate.  The  ratio  r)rm,/rirms  in  (48) 


scales  the  energy  in  the  time  series  r(n  to  match  that  of  ijm . 

The  directional  frequency  spectra  of  the  77"  time  series  are  computed  for  values  of  the  constant  c  of 
0.33  and  1.0,  with  the  results  shown  in  Figures  11  and  12.  Additionally,  Figure  9  may  be  interpreted  as  the 
c  =  0  case.  As  c  increases,  the  angular  domains  of  the  0.3  and  0.1  m2  Hz”1  rad"1  contours  increase  in  the 
clockwise  direction.  Asymmetry  is  evident  even  at  the  c  =  0.33  level.  Figure  12  should  correspond  to  the 
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sum  of  Figure  9  and  Figure  10  (after  rescaling)  rotated  by  120*  if  the  time  series  rj  and  rf  are  uncorrelated, 
which  is  true  in  this  case.  The  true  peaks  due  to  r,  and  if  should  be  located  at  about  145°  and  280°, 
respectively,  while  Figure  12  shows  two  peaks  at  about  160°  and  255°.  Despite  the  slightly  decreased 
angular  separation  between  the  two  peaks  in  the  latter  case,  the  presence  of  two  peaks  in  the  c  =  1  case 
shows  that  the  generalized  DASE  method  is  capable  of  resolving  a  bimodal  spectrum  with  similar  peak 
energy  and  peak  angles  separated  by  125°.  Furthermore,  the  asymmetry  of  contours  in  the  c  =  0.33  cases 
indicates  that  the  method  is  sensitive  enough  to  respond  to  the  presence  of  significantly  smaller  secondary 
waves  separated  in  angle  from  dominant  waves.  This  allows  the  generalized  DASE  method  to  be  used  in 
detecting  the  energy  of  reflected  waves  from  the  RASEX  tower  posts  (Hare,  1995). 
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6  Application  of  Motion  Correction  Scheme  and  Extended  DASE 
Method 

In  this  section  we  combine  the  motion  correction  scheme  discussed  in  Section  2  and  the  extended  DASE 
method  developed  in  Section  4.  Using  signals  obtained  from  a  towed  research  catamaran  during  CoOP 
experiment,  directional  wave  spectra  are  estimated. 

v 

6.1  Overview  of  CoOP  experiment 

The  Coastal  Ocean  Processes  (CoOP)  experiment  occurred  in  April/May  of  1995,  off  the  California  coast, 
and  was  sponsored  by  the  National  Science  Foundation.  The  principle  objective  of  the  program  was  to 
investigate  the  role  of  various  ocean  surface  processes  on  the  air-sea  gas  exchange.  The  research  catamaran 
was  towed  to  the  side  of  the  R/V  New  Horizon.  Measurements  of  gravity-capillary  waves,  surface  gravity 
waves,  sub-surface  currents/turbulence,  and  surface  enrichment  were  made  possible  with  the  instrument 
array  aboard  the  catamaran.  This  array  included  a  set  of  six  capacitive  wave  wires,  a  motion  detection 
package,  a  three-axis  acoustic  Doppler  current  meter,  a  hot  film  anemometer,  and  a  surface  skimmer.  The 
motion  detection  package  consisted  of  a  three-axis  translational  accelerometer  and  a  three-axis  angular- 
rate  sensor.  Both  the  wave  wire  data  and  the  motion  detection  package  data  were  recorded  digitally  at  a 
sampling  rate  of  100  Hz.  The  wind  speed/direction,  and  the  wind  stress  were  obtained  at  the  bow  of  the 

ship  using  a  sonic  anemometer. 

6.2  Directional  frequency  spectra  results 

The  procedure  for  estimating  the  wave-height  spectrum  for  CoOP  was  identical  to  the  procedure  used 
during  High- Res  2.  Because  no  independent  measurements  of  the  wave-height  spectrum  were  available,  we 
cannot  validate  the  spectra  for  CoOP  as  we  did  for  High-Res  2.  In  Figures  13  through  15,  the  results  of 
frequency  spectra,  directional  frequency  spectra,  and  normalized  directional  spreading  are  shown  for  three 
cases.  Wave-height  frequency  spectra  for  Figures  13-15  were  computed  using  102  s  time  blocks  and  total 
time  series  of  39,  46,  and  45  minutes,  respectively.  All  measurements  were  performed  in  the  vicinity  of 
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Santa  Catalina  Island,  shown  in  Figures  16.  The  first  two  figures  correspond  to  cases  where  wind  and  swells 

are  aligned.  Figure  13  shows  an  example  of  well-developed  swells  and  relatively  low  local  wind  waves  under 

weak  wind.  The  directional  spreading  of  swells  (0.1-0.2  Hz)  is  narrower  than  wind  waves  (0.3-0.4  Hz)  as 

expected.  On  the  other  hand,  the  wave  field  in  Figure  14  is  mostly  local  wind  waves  generated  by  strong 

wind.  The  low  level  of  swells  may  be  partially  due  to  the  sheltering  effect  of  the  island.  In  both  figures, 

the  directions  of  local  wind  waves  are  consistent  with  the  observed  wind  directions. 

\ 

An  example  is  shown,  in  Figures  15,  where  swells  and  local  wind  waves  are  not  aligned.  In  this 
case,  the  wind  speed  is  relatively  low  and  both  wind  speed  and  wind  direction  are  rather  variable.  At  swell 
frequencies  (0. 1-0.2  Hz)  the  directional  spreading  is  narrow  and  well-defined.  At  the  highest  frequency  (0.4 
Hz),  the  peak  direction  is  roughly  consistent  with  the  local  wind  but  the  wave  energy  is  more  widely  spread 
in  all  directions.  In  between  at  0.3  Hz,  the  wave  field  appears  to  be  confused  because  of  the  presence  of 
both  swells  and  local  wind  waves. 

These  examples  demonstrate  that  the  combination  of  our  motion  correction  scheme  and  the  ex¬ 
tended  DASE  method  may  successfully  estimate  directional  wave  spectra,  including  both  swells  and  local 
wind  waves  in  the  range  of  wave  conditions  encountered  during  CoOP.  A  more  comprehensive  analysis  of 
the  accuracy  of  this  method  in  a  wider  range  of  conditions  is  a  topic  for  future  research. 
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7  Conclusion 


We  have  developed  a  method  to  obtain  directional  spectra  of  surface  gravity  waves  using  a  motion  detection 
package  and  a  wave  wire  array,  both  mounted  on  a  moving  platform.  The  validity  of  the  motion  correction 
scheme  to  the  wave  wire  signals  has  been  confirmed  by  comparison  with  independent  buoy  observations. 
An  improved  method  of  computing  the  directional  frequency  spectrum  has  been  derived  by  generalizing 
the  DASE  method  of  Davis  and  R^gier  (1977)  to  allow  for  the  motion  of  the  observer  relative  to  the  wave 
field.  The  applications  of  the  method  to  two  field  data  sets,  one  from  a  fixed  platform,  and  the  other  from 
a  moving  platform,  have  shown  that  the  propagation  direction  of  the  dominant  wind  waves  agree  well  with 
the  wind  direction.  Furthermore,  the  swell  direction  at  lower  frequencies  can  be  resolved  accurately  even 
if  it  is  not  aligned  with  local  wind.  A  numerical  test  with  the  fixed  platform  data  demonstrates  sufficient 
sensitivity  of  the  model  to  detect  the  presence  of  secondary  waves  propagating  at  different  incidence  angles 
than  those  of  dominant  waves. 

The  directional  wave  information  obtained  in  this  study  can  be  directly  correlated  with  other 
collocated  physical  measurements  of  atmospheric  turbulence,  short  wind  waves,  and  subsurface  turbulence 
in  order  to  address  detailed  interaction  between  various  air-sea  interaction  processes  and  surface  gravity 
waves.  Currently,  we  are  using  the  directional  spectra  of  wind  waves  in  a  study  of  their  modulation  of 
shorter  gravity-capillary  waves. 
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Appendix  A:  Estimation  of  Frequency  Spectra 

The  Welch  method  is  a  standard  means  of  computing  cross-spectra,  and  utilizes  overlapping  time  blocks 
and  non-rectangular  data  windows  (Marple,  1987).  The  Welch  estimate  of  the  cross-spectrum  is  given  by 

Cn(fk)  =-±p'£Xi(fk,P)X;(fk,P)  .  (49) 

where  *  denotes  the  complex  conjugate  operation,  the  summation  is  over  P  time  blocks,  each  of  duration 
T,  and  Xx  and  X2  are  the  Fourier  transforms  of  the  time  signals  xx(t)  and  x2(t).  If  x2(t)  =  xx(t)  then 
the  result  Cxx(fk)  is  the  autospectrum  of  xx(t).  The  constant  p  is  described  shortly.  The  discrete,  finite 
Fourier  transforms  (DFT’s)  Xx{fk,p)  and  X2{fk>p)  are  given  by  the  general  relation 

X(A)  =  A  E  .  m 

Jsamp  n=0 

where  N  is  the  number  of  discrete  points  in  the  time  block,  f,amP  is  the  sampling  frequency  in  Hz,  and 
x(tn)  is  the  discrete  time  series.  The  discrete  frequencies  fk  are  given  by 

fk  =  f,amPk/N  with  k  =  0,1,...  ,JV- 1  •  (51) 

In  practice,  the  DFT’s  are  computed  using  the  fast  Fourier  algorithm  (see  e.g.,  Bendat  and  Piersol,  1986  ). 
We  have  generally  chosen  w(tn)  to  be  the  Hann  window  (Marple,  1987). 

,w-«*(^^!ia).  («> 

This  window  has  smaller  sidelobes  than  those  of  the  simple  rectangular  window  u»r(t„)  =  1.  Sidelobes  are 
a  phenomenon  of  discrete  Fourier  analysis  causing  the  power  at  one  frequency  to  be  spread  to  neighboring 
frequencies.  Although  the  Hann  window  has  smaller  side  lobes  than  the  rectangular  window,  it  has  a 
slightly  wider  mainlobe  peak.  This  mainlobe  width  determines  the  effective  resolution  of  the  spectrum. 
The  rectangular  window  has  an  effective  mainlobe  width  equal  to  the  frequency  spacing  of  (51),  or  Sf  = 
f\  =  f samp IN.  For  the  Hann  window,  this  resolution  is  1.5  5f. 

The  constant  p  in  (49)  compensates  for  the  loss  in  total  energy  of  the  time  signal  a;(t„)  due  to  the 

windowing  operation.  It  is  defined  for  all  windows  by 

„  =  i  £«,(<„).  M 

71  =  0 
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For  the  Hann  window,  fi  =  0.37. 

The  use  of  overlapping  time  blocks  in  the  calculation  of  Cnifk)  increases  P  and  therefore  the 
statistical  stability  of  the  estimate.  It  is  justified  by  the  form  of  the  Hann  window.  Consider  two  successive 
time  blocks  of  length  N.  Signals  near  their  common  edge  are  damped  heavily  by  Hann  window  values  close 
to  zero,  as  shown  by  (52),  with  n  approaching  0  or  N  -  1.  Thus  their  contribution  to  the  DFT  of  (50)  is 
small.  In  an  overlapping  block  centered  on  this  edge,  these  oscillations  are  windowed  with  wh{tn)  close  to 

N 

1,  and  will  make  a  significant  contribution  to  the  DFT.  The  sum  of  the  window  values  u>h(t„)  of  one  block 
and  the  values  of  the  2  overlapping  windows  is  in  fact  identically  equal  to  1.  Thus,  only  the  first  and  last 
N/2  points  of  the  entire  time  series  are  underutilized. 
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Figure  Captions 

Figure  1:  Schematic  showing  wave  height  estimation  from  the  research  catamaran. 

Figure  2:  Flow  chart  of  the  wave-height  estimation  procedure. 

Figure  3:  The  location  of  the  High-Res  2  experiment,  near  the  barrier  islands  of  North  Carolina.  Ship 
locations  are  shown  as  crosses.  The  East  and  West  discus  buoys  are  indicated  by  filled  circles.  The  average 
location  of  the  Gulf  Stream  North  wall  is  shown. 

Figure  4:  High-Res  2  wave-height  frequency  spectrum  for  Julian  day  171  from  15:10  to  17:11  GMT  ( P  =  60 
time  blocks),  u.  =  0.21  m  s"1  and  U10  =  6.4  m  s'1.  Solid:  without  Doppler  shifting  correction.  Dotted: 
with  unidirectional  Doppler  shifting  correction.  Dashed:  West  discus  buoy  spectrum  at  16:00  GMT. 

Figure  5:  High-Res  2  wave-height  frequency  spectrum  for  Julian  day  173  from  19:18  to  20:48  GMT  (P  =  45 
time  blocks),  u.  =  0.24  m  s"1  and  Uio  =  8.0  m  s"1.  Solid:  without  Doppler  shifting  correction.  Dotted: 
with  unidirectional  Doppler  shifting  correction.  Dashed:  East  discus  buoy  spectrum  at  20:00  GMT. 

Figure  6:  Solutions  to  the  Doppler-shifted  dispersion  relation  (13)  for  deep  water  and  U0bs  =  0.8  m  s1. 
Dotted  lines  show  the  case  of  /  =  ±0.2  Hz.  The  k  solutions  for  /  =  ±0.2  Hz  are  given  by  equations  (14) 

and  (15). 

Figure  7:  ka,  kp,  and  klt  solutions  to  the  Doppler-shifted  dispersion  relation  (13)  for  deep  water  and 
fixed  /  of  0.29,  0.38,  0.48,  and  0.58  Hz.  Wavenumber  magnitude  k  increases  in  the  radial  direction  and 
propagation  direction  9k  increases  in  the  counterclockwise  direction  from  the  positive  x  axis.  U0bs  =  0.8  m 
s-1,  dobs  =  0°.  Together  the  three  solutions  form  a  closed  loop  only  for  /  <  /.  =  0.49  Hz. 

Figure  8:  The  location  of  the  RASEX  experiment,  near  the  coast  of  Denmark.  The  research  tower  shown 
was  outfitted  with  an  array  of  three  wave  wires. 

Figure  9:  RASEX  wave-height  spectral  estimates  for  Julian  day  121  from  5:52  to  6:22  GMT  (P  =  34  time 
blocks),  (a)  Frequency  spectrum  in  m2  Hz'1,  (b)  Directional  frequency  spectrum,  with  intrinsic  frequency 
(scale  shown  on  axes)  in  the  radial  direction  and  wave  propagation  angle  (North  upwards,  East  to  the  right) 
in  the  tangential  direction.  Contour  labels  are  in  units  of  m2  Hz-1  rad-1  and  adjacent  contours  have  ratio 
101/2.  Arrow  represents  wind  direction;  wind  speed  is  12  m  s-1. 

Figure  10:  RASEX  wave-height  spectral  estimates,  as  in  Figure  9,  but  for  Julian  day  121  from  10:35  to 
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11:05  GMT  ( P  =  34  time  blocks).  Wind  speed  is8ms". 

Figure  11:  Test  of  the  angular  resolution  of  the  DASE  method  (P  =  34  time  blocks).  Spectral  estimates  of 
rf*  in  equation  (48),  with  c  =  0.33.  Parts  (a)  and  (b)  as  in  Figure  9. 

Figure  12:  Test  of  the  angular  resolution  of  the  DASE  method  ( P  =  34  time  blocks).  Spectral  estimates  of 
r f  in  equation  (48),  with  c  =  1.0.  Parts  (a)  and  (b)  as  in  Figure  9. 

Figure  13:  CoOP  wave-height  spectral  estimates  for  Julian  day  132  from  15:46  to  16:25  GMT  (P  —  43 

V. 

time  blocks),  (a)  Frequency  spectrum  in  m2  Hz-1.  Solid  line:  intrinsic  frequency,  dotted  line:  observed 
frequency,  (b)  Directional  frequency  spectrum,  with  intrinsic  frequency  (scale  shown  on  axes)  in  the 
radial  direction  and  wave  propagation  angle  (North  upwards,  East  to  the  right)  in  the  tangential  direction. 
Contour  labels  are  in  units  of  m2  Hz-1  rad-1  and  adjacent  contours  have  ratio  101/2.  Arrow  represents 
wind  direction;  wind  speed  is  4.5  m  s"1.  (c)  Normalized  directional  spreading  at  different  frequencies.  Solid 
line:  0.1  Hz;  dotted'  line:  0.2  Hz;  dashed  line:  0.3  Hz;  dash-dot  line:  0.4  Hz. 

Figure  14:  CoOP  wave-height  spectral  estimates  for  Julian  day  133  from  18:44  to  19:30  GMT  (P  =  52 
time  blocks),  (a)  Frequency  spectrum  in  m2  Hz-1.  Solid  line:  intrinsic  frequency,  dotted  line:  observed 
frequency,  (b)  Directional  frequency  spectrum,  with  intrinsic  frequency  (scale  shown  on  axes)  in  the 
radial  direction  and  wave  propagation  angle  (North  upwards,  East  to  the  right)  in  the  tangential  direction. 
Contour  labels  are  in  units  of  m2  Hz”1  rad-1  and  adjacent  contours  have  ratio  101/2.  Arrow  represents 
wind  direction;  wind  speed  is  9  m  s-1.  (c)  Normalized  directional  spreading  at  different  frequencies.  Solid 
line:  0.1  Hz;  dotted  line:  0.2  Hz;  dashed  line:  0.3  Hz;  dash-dot  line:  0.4  Hz. 

Figure  15:  CoOP  wave-height  spectral  estimates  for  Julian  day  131  from  04:09  to  04:54  GMT  ( P  =  51 
time  blocks),  (a)  Frequency  spectrum  in  m2  Hz-1.  Solid  line:  intrinsic  frequency,  dotted  line:  observed 
frequency,  (b)  Directional  frequency  spectrum,  with  intrinsic  frequency  (scale  shown  on  axes)  in  the 
radial  direction  and  wave  propagation  angle  (North  upwards,  East  to  the  right)  in  the  tangential  direction. 
Contour  labels  are  in  units  of  m2  Hz-1  rad-1  and  adjacent  contours  have  ratio  101/2.  Arrow  represents 
wind  direction;  wind  speed  is  4  m  s"1.  (c)  Normalized  directional  spreading  at  different  frequencies.  Solid 
line:  0.1  Hz;  dotted  line:  0.2  Hz;  dashed  line:  0.3  Hz;  dash-dot  line:  0.4  Hz. 

Figure  16:  The  location  of  the  CoOP  wave  measurements  of  this  study,  near  Santa  Catalina  Island. 
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